% Calculate input impedance 

aantal=length(pressure);
T=max(time);

% compute input impedance

Zin=fft(pressure)./fft(flow);
Zinmod=abs(Zin);
Zinphase=unwrap(angle(Zin));
freq=(1/T)*((1:aantal)-1);

figure;
subplot(2,1,1);
title('modulus input impedance');
xlabel('frequency (Hz)');
ylabel('impedance (mmHg/(ml/s))');
semilogy(freq(1:11),abs(Zin(1:11)),'m+',freq(1:11),abs(Zin(1:11)),'b');
subplot(2,1,2);
title('phase input impedance');
xlabel('frequency (Hz)');
ylabel('phase angle (°)');
plot(freq(1:11),180*angle((Zin(1:11)))/pi,'m+',freq(1:11),180*angle((Zin(1:11)))/pi,'b');

fprintf(1,'Harmonic\tFrequency\tModulus\tPhase\n');
for k=1:11
   fprintf(1,'%g\t%g\t%g\t%g\n',k,freq(k),abs(Zin(k)),(180/pi)*angle(Zin(k)));
end   
